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A quasi 2-dimensional recursive lattice formed by planar elements have been designed to investi¬ 
gate the surface thermodynamics of Ising spin glass system with the aim to study the metastability 
of supercooled liquids and the ideal glass transition. The lattice is constructed as a hybrid of 
partial Husimi lattice representing the bulk and ID single bonds representing the surface. The 
recursive properties of the lattices were adopted to achieve exact calculations. The model has an 
anti-ferromagnetic interaction to give rise to an ordered phase identified as crystal, and a metastable 
solution representing the amorphous/metastable phase. Interactions between particles farther away 
than the nearest neighbor distance are taken into consideration. Free energy and entropy of the 
ideal crystal and supercooled liquid state of the model on the surface are calculated by the partial 
partition function. By analyzing the free energies and entropies of the crystal and supercooled liquid 
state, we are able to identify the melting transition and the second order ideal glass transition on the 
surface. The results show that due to the coordination number change, the transition temperature 
on the surface decreases significantly compared to the bulk system. Our calculation agrees with 
experimental and simulation results on the thermodynamics of surfaces and thin films conducted 
by others. 


I. INTRODUCTION 

Glass transition on surface/interface/thin film has drawn intensive interests in the last two decades for two reasons. 
Firstly, the importance of surface and thin film in materials science and engineering requires a better understanding 
of its dynamic and thermodynamic properties. Secondly, the confined geometry is a good approach to understand 
the mysterious glass transition itself, especially to explore the dynamic properties within the thin film geometry. 
Numerous works of both experimental and simulation/calculation approaches have been done [1-23] on glass transition 
on surface/thin film. 

Keddie and co-workers firstly investigated the supported thin film of PS by ellipsometric measurements Q- They 
prepared several PS films supported by silicon wafers with thickness from 100 A to 3000 A. The measurements 
indicated the deduction of the Tg with thickness under 400A. An empirical relationship between thickness and Tg 
was given as: 


Tg = rB-"^[i-(^)^], (1) 

where is the glass transition temperature of the bulk PS. The a and S are the parameters fitted to be 32 A and 

1.8 respectively, h is the thickness of the film. Following Keddie’s work, many researches have been done on different 
supported thin films of various polymers [2-3] by different characterization methods, such as X-ray reflectivity, positron 
annihilation, and dielectric [^-Ig . Most results have demonstrated the same phenomenon that for liner polymer the Tg 
decreases with the thickness of films. However the supported film has a considerable film-substrate interaction, which 
makes the conclusion controversial. Strong attractive interaction between the substrate and thin film may increase 
the Tg of thin film above the bulk Tg [6]. van Zanten et al. measured Tg of poly-2-vinylpyridine on oxide-coated 
Si substrates, and found it increase by 50K than the bulk Tg, for a 77 A film [7]. Forrest and co-workers have done 
pioneer works in measuring the Tg of free-standing thin films [8-10]. They measured the Tg of free-standing PS films 
with thickness from 200 to 2000 A and different molecular weights by Brillouin Light Scattering and transmission 
ellipsometry. Their results showed that the Tg decreases with the thickness of PS thin film with a much larger 
magnitude: for example, the Tg of 200 A film with Mw within 120K to 378K reduces by 70K below the while 
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FIG. 1: A regular square lattice and its analog of a recursive lattice (Husimi lattice). 


this magnitude is around lOK for supported films. The empirical equation (EqlT]) derived for supported film still holds 
for the low Mw free-standing films. With 6 = 1.8, the parameter a was found to be 78 A which is twice of it found 
for supported films. 

Other than the experimental work, computer simulation and calculation have also been developed to investigate the 
glass transition on surface/thin film, and most of them are for polymer systems [11-23]. Molecular Dynamics (MD) 
and Monte Carlo (MC) method were usually employed with various modelings, and confirmed the Tg decrease with 
the thickness reduction for both supported and free standing film, or Tg increase in some particular substrate-film 
cases. The EqU) derived from experiments can also be validated by simulations, nevertheless the explanation for the 
mechanism of Tg reduction is still a matter of debate. Most MD simulations verified the experimental observation 
that for supported film, a strong substrate-film interaction will increase the Tg above the value in the bulk, while the 
weak substrate-film interaction will lead Tg reduction, and free-standing film shows a much larger Tg reduction than 
supported film with weak substrate-film interaction [12]. Mattice and co-workers firstly reported the MC simulation in 
this field [13, 14]. de Pablo et.al. reported a MC simulation on free-standing films of both linear and cyclic polymeric 
chains [23]. Basically the Tg reduction with smaller thickness was confirmed in the above investigations. In de Pablo 
and co-workers’ work, the relation between Tg and can also be fitted to be Eq. [T]with slight difference on the 

parameters fitting. 

Although the glass transition is not a unique property of polymers, most simulation works were modeled for polymer 
systems, while the works on small molecule systems are very rare. Meanwhile, Ising spins glass has also been widely 
utilized to study the glass transition [24-37]. By different lattices adopted and interactions setup the Ising model is 
capable to describe various systems, such as gas, liquid, crystal and glass, and consequently the phase transitions, 
like melting and glass transition [24]. However, very few of the efforts were on surface/thin film glass transition of 
Ising spin glass. This field had been explored by Gujrati et al. by applying Ising model on a modified Bethe lattice 
to describe the thermodynamics of polymer systems near surface [38-42]. In this paper we follow the similar method 
to study the glass transition of Ising spin glass on the surface/thin film on a specially constructed recursive lattice. 


II. SURFACE RECURSIVE LATTICE (SRL) GEOMETRY 

Except in some rare cases [43-45], a many-body system with interactions on a regular lattice is difficult to be solved 
exactly because of the complexity involved with treating the combinatorics generated by the interaction terms in the 
Hamiltonian when summing over all states. The mean-field approximation is commonly adopted to solve this problem, 
e.g. the Flory model of semiflexible polymers [46, 47]. On the other hand, recursive lattices enable us to take the 
explicit treatment of combinatorics on these lattices and no approximation is necessary [25-28, 48]. The recursive 
lattice is chosen to have the same coordination number as the regular lattice it is designed to describe. As usual, 
the coordination number q is the number of nearest-neighbor sites of a site. A typical recursive lattices, the Husimi 
lattice, as the analogs of the 2-D lattice is shown in Fig. [TJ For the bulk system, the recursive lattice calculations 
have been demonstrated to be highly reliable approximations to regular lattices [48-50]. 

In this work we construct a recursive lattice to describe the surface/thin film. This surface recusive lattice (SRL) 
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FIG. 2: A regular 2D square lattice of a thin film with thickness = 5. 


is to mimic the 2D case, i.e. the 2D bulk with ID surfaces. The SRL is integrated of square units representing the 
bulk and single bond representing the surface, the structure is made to have the same coordination numbers with 
regular 2D square lattice. Ising spins will be applied on the lattices to represent a small molecule system. The exact 
calculation will be drived and the solutions will be discussed. 


A. Construction of SRL 

1. Structure 

To approximate the regular lattice due to the same coordination number, we firstly need to figure out the co¬ 
ordination number and interactions of a regular lattice of surface/thin film to construct the SRL to represent the 
surface/thin film. The Fig. [5] shows a regular 2D square lattice of a thin film with thickness equals to 5 and the thick 
line is the surface. (In all figures of this paper the thick bond represents surface bond.) Here we label the central 
layer to be the 0th level, and the layer next to 0th layer is the level 1st. The surface layer is labeled as level S. The 
sites on the surface have coordination number of 3, while inside the bulk the coordination number is 4. That is, a 
hybrid recursive lattice with q = 3 and 4 is required to describe the surface/thin film. 

Now we take out the basic units inside the bulk and on the surface (Fig. [3]) to construct a recursive lattice. Inside 
the bulk we simply have Husimi lattice of g = 4. For the surface structure, a single bond representing in the surface 
unit links on a square unit then we can have sites with q = 3. However this unit cannot be simply adopted to 
construct a recursive structure, because the recursive calculation technique (will be discussed later) requires an origin 
point to which the entire tree is symmetrical. For the unit shown in Fig. 0^, wherever we determine the origin is 
(point A or B), the unit is not symmetrical to that point. We modify the surface unit by replacing one square by an 
artificial single bond at one lower corner. Although this unit differs from the regular lattice we want to approximate, 
the coordination numbers on the surface square still accord to our design and the calculation in the following sections 
will show that this approximation is practical. The modification is shown in Fig. 01 the modified unit (Fig. Hb) is 
then symmetrical to the point A. 
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FIG. 3: The basic unit in the bnlk and on the surface of a regular 2D square lattice 
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FIG. 4: The modification of surface square unit. The bulk square on lower right corner is replaced by a single bond. In this 
way the modified structure is symmetrical to an imaginary axis links point A and the center of square unit S. 


Then we can put the bulk and surface units together to construct a recursive lattice to approximate the regular 
thin film lattice. The structure of a thickness = 5 lattice is shown in Fig. [5] 

From surface A to B, we have a finite bulk portion with a thickness = 5, and these five layers are labeled as S', 1, 
and 0 like the regular film lattice in Fig. [5] The thick single bond representing surface links two identical bulk trees. 
In this example, one bulk tree Tb is linked overall with another 36 trees (only one bulk tree is drawn in above Fig. 
ED, while they are all identical and independent with each other except a single surface bond connection. Recursively, 
each of these 36 trees is also linked with another 35 trees by single surface bonds. This lattice is an infinite tree 
integrated by the finite size bulk portions and infinitely large surfaces. If we look at the surface integrated by thick 
surface bonds indicated by the arrow at the right lower corner in Fig. [51 we can see the surface going through the 
bulk trees Tb is infinitely large with a coordination number 3 everywhere. A stretched surface is shown in Fig. [5] to 
provide a more obvious view of the surface. Comparing to the regular lattice surface, the main difference is that one 
surface also receives thermodynamic contributions from other surface structures, which is caused by the surface unit 
modification. But this is an approximation we have to take to achieve the recursive calculation. Although we have 
infinite number of bulk trees and surfaces in this lattice, since they are independent and identical, it does not impact 
the thermodynamic properties of a local region we are going to investigate. 
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2. The sites labeling on bulk and surface units 
We label each site in a square as shown in Fig. [7j 

The site in a square unit is labeled as Sa,b, where a is the index of the square unit in the bulk tree, b is the index 
of site. Therefore, one site Sa,b actually has two labels depending on which square unit we refer to. For example the 
site Sa ,4 is also S’(o+i)_i in the higher level unit. When we are only focusing on one unit, the site Sa,b can be denoted 
as Sb for convenience. In a square unit, the base site is determined to be the site closest to the bulk tree origin (The 
0th square) and labeled as Sa,i- Therefore, in the upper enlarged circle in Fig. [3 the base site is the lowest site, while 
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FIG. 6: The stretched view and the contribution direction of one surface in SRL. 


in the lower enlarged circle, the base site is the left one. The labeling in the origin unit is a special case which will 
be discussed later. The interactions are also sampled in the enlarged circle: J is the interaction between the nearest 
sites, Jp is the interaction between the second-nearest sites, i.e. the diagonal sites. 

For surface calculation, we label the surface sites as shown in Fig. [8l Due to our calculation method, only two 
labels are necessary in surface labeling. The site close to the bulk site is labeled as S'o, and the site being diagonal to 
the bulk site is labeled as Si. 

1.3) The interactions and energy in bulk square unit 

We consider four kinds of interactions in our model: J the interaction energy between the nearest sites, Jp the 
interaction energy between the second nearest sites, J/ the interaction energy of three sites (triplet), and J” the 
interaction energy of four sites (quadruplet). 

We introduce the following variables to count the interactions and magnetic fields of one square unit: 


klnb = Si X 5*2 “t“ *5'i X Sn -f S 2 X -p 5*3 X 5*4; 

Ad = Si X Si -\- S2 S^\ 

Atri ~ y ‘5*2 ^ *^3 T ^ ^2 ^ ^4 T X S^ X Si -f S 2 X S^ X 5*4; 

Aqua y ^2 y *^3 ^ *^4: 

Amag = S2 + S3 + Si- 

Where Sm is an Ising spin and has the value of +1 or —1. Note that the base site is not included in the magnetic 
term Amag, because the calculation has a recursive fashion that the energy of one unit is included in the contribution 
to the next level’s unit, and the base site of one unit will be counted as the top site in the next level’s energy equation, 
thus we need to exclude it here to avoid the double counting. 

Then, for a particular cell a with a certain spin conformation F, the energy of the Ising model in one cell is: 


e„ = - J • Aub(r) - Jp • Ad(r) - J/ • Atri(r) - j” • Aqua(r) - ff ■ Amag(r). 


( 2 ) 


The total energy of the Ising model on SRL is the sum of the energy of all cells: 
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Etotai = ^e„(r). (3) 

OL 

In this work we set J to be negative to adopt the the anti-ferromagnetic Ising model, the neighbor spins preferring 
different states corresponds to the lowest energy stable state, i.e. the ideal ordered crystal conformation, which have 
the largest weight comparing to other conformations at the same temperature, while the neighbor spins of the same 
states are the most unstable conformations with the highest energy, which represent the disordered phase. 

Then the Boltzmann weight u;(r) of state T is given by 


cc(r)=exp[-^(e(r)-i7,5o)], (4) 

where /? = l/feeT is the inverse temperature, and here we set the Boltzmann constant fee = 1 to make the 
temperature to be generalized in the unit of energy. 

Then the partition function of a finite lattice with a cells is the sum over products over the Boltzmann weight: 

Z{T)= Y. [n^«(r)]. 

{S=±l} a 


( 5 ) 
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FIG. 8: The surface sites labeling in SRL. 


3. The interactions and energy on surface 

In Fig. [BJ the single bond unit and the square unit alternatively appear on the surface structure. The interactions 
and energy of the surface square unit is similar to the bulk square unit as we discussed in the previous section. 
However, here we would like to assign a different nearest neighbor interaction J on the surface bond, to distinguish it 
with J in the bulk. This enables us to investigate more complex surface properties. And similarly we may also have 
a different diagonal interaction (Jp), triplet interaction (J/), quadruplet interaction (J”) and magnetic field (H): 


e„ = - J . {S2 XS4 + S3XS4)-J- {Si xS2 +Six S3) - Jp ■ Ad(r) - Jt ■ Atri(r) - J” • Aqua(r) - H ■ A„,ag(r). (6) 
On the surface single bond unit, the interaction is much simpler since there is only one nearest neighbor interaction: 


e=-J-{SoxSi)-H-Si. (7) 

Similar to the role of base site in bulk calculation, here the magnetic field of So is not included in the second term 
because that site will be counted in the next level’s energy equation. The ‘level’, in this statement, is indexed by 
taking the direction indicated by the arrow in Fig. |6l and employing an imaginary origin point infinitely far away 
from the region we concern. Since the surface is infinitely large, the selection of ‘origin point’ does not affect us to 
achieve the solutions on surface. 

By the setup of interaction energy parameters, we can simulate various systems with particular interactions and 
energy to study their thermodynamic properties and phase transition with the determination of the partition function. 
The effect of energy parameters will be discussed in section IV. We will generally set J = —1 to determine the 
temperature scale for our antiferromagnetic model. The solution based on J = — 1 and all other parameters to be 0 
is called the reference model. 


B. General recursive calculation technique 

To discuss the calculation of our model, we first need to introduce the concept of sub-tree contributions. In the 
finite bulk tree we have an origin at the center of the tree. For each square unit, the base site is the closest site to the 
origin, and there are three sub-trees coming from the other three sites. The sub-trees could either be three identical 
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portions of the bulk tree, or three surface trees linked by the single surface bond if it is the square unit on the surface. 
We label the levels in one unit and show the sub-tree contributions in Fig. 

By this index, the sites on the level m is represented as Sm- The sub-tree going to the site on level m is marked as 
Tm- In this way, we can introduce the partial partition function (PPF) Zm^.i(Sm+i) at the level m -|- 1 to represent 
the contribution of the branch Tm+i with a certain state of spin Sm+i as its base site to the lower level m, and Tm 
to the lower level m — 1, etc. Therefore the PPF Zm{Sm) at level m is a function of the PPFs Z^+iiSm+i) and 
Zm+2{Sm+2) at higher level m -I- 1 and 771 + 2 and the local weight. Then we can start from the highest level, count 
the contributions of sub-trees on each level and recursively go to the lower level unit the origin point, where it gains 
the contribution of the entire lattice, and the thermodynamics of the entire system can be obtained by the partition 
function. The PPF of a sub-tree Tm ^m(±) is the sum of the configurations with the spin Sm = ±1 over the products 
of the PPFs on higher levels with the local weight a;(r). For a square unit, Zm{+) has 8 terms with Sm = +1 and 
Zmi—) is the sum of the other 8 configurations with Sm = — 1: 

8 2 

Zn.{+) = Y.^ n Zm+l{.Sm+l)]Zm+2{Sm+2)w{Ts^ = l), (8) 

F—1 {m+1} 

8 2 

^™(-)=E[n ^m+l(<S'm+l)]^m-|-2(<5'm+2)u;(F5^=_l), (9) 

F—9 {m+1} 

Depending on the spin state (S'o = ±1) of the origin site, the two sub-trees’ contributions to the whole system are: 

^o(+)exp(/3ff) or Z^(-) exp(-^iJ). 

where the magnetic term is to count the contribution of origin site which is not contained in Eql^ The partition 
function Zq of the entire lattice can be accounted by the contributions to the origin site: 


Zo = Z^(-h)exp(^iJ)-hZ^(-)exp(-/3+. (10) 

The first term is to count the conformation with the origin site spin Sq = -1-1, the square of PPF represents two 
sub-trees and the weight of the origin spin itself is counted as exp(/liJ), similiarly for the S'o = — 1 situation in the 
second term. 

Or, if the origin is defined to be a square unit, the partition function is 

16 4 

Z77=Y.{X{Z,{S,)-U7{T)) 

F i 
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where Si is one of the four spins on the origin square, Zi[Si) is the PPF of the sub-tree coming into the site 5”^, 
and w(r) is the local weight of the origin square. 

Then we can introduce the ratios 


— 


Zm{+) 


^m( + ) + Zm{ — ) 


t Um — 


ZU-) 


Zm{S-) + Zm( — ) 


( 11 ) 


These two are the ratios of PPFs on the level m; they indicate, although not exactly to be, the probability that 
the site Sm on level m is occupied by a plus spin (xm) or a minus spin (j/m)- Therefore, these two ratios, which we 
call the solutions of the lattice, are critical in our model and all the thermodynamic calculations are based on these 
solutions. 

Define a compact notation for convenience 


^m{Sm.) — 


if Sm — +1 

if Sm = -1 


( 12 ) 


by introducing 

Bm = Zm{-\-) + Zm{ — ), 

we then have Zmi+) = BmXm and Zm(—) = BmUm- Substituting these two equations into Eql8] and |9] gives 


ZmiSm) = Zm+l{Sm+l)]Zm+2{Sm+2)w{T), 

r {m+1} 

where the sum is over F = 1, 2,3,..., 8 for Sm = +1, and over F = 9,10,11,..., 16 for Sm = —1- It is clear that 
the solution (ratios) on one level is a function of the solutions on higher levels in a recursive fashion. 

For convenience, we introduce the polynomials: 


8 2 

Etn (‘5'm+l)]'2(m+2 ('5^m+2 )'^(r) j 

P—1 {m+1} 

8 2 

Qm—(ym+lj ym+ 2 ) — Etn (‘5'm+l)] Zm+2{Sm+2)u){T), (14) 

P—9 {m+1} 

Qm (^m+1 1 ^m+2) — Qm-\- (-^m+l: •^m+2 Qm— (^m+l j ym+ 2 ) — ~^2 d ' 

-^m+l-^^n+2 

Then the ratio Xm is just as function of the ratios on higher level x^+i and Xm-\- 2 ' 


Q-\-(^Xm-\-l': Xm-\-2') 
Q(^m+1; ^m+ 2 ) 


(16) 


1. Fix-point solution 

From the recursive relation shown in Ea llhl we may expect a repeating solution, which has the form: 
xl, x2,..., Xy{Xv > 1) and Xk = Xy+k for some value of u > 1, 
so that Eg fT6l holds. 

Such a solution can be called a ?;-cycle solution that repeats after the w-th application of the recursive relation 
equation. 

For example, for a 1-cycle solution x, we simply have the same ratio x on all the levels 


Q+{x,x) 

Q{x,x) 


X = 


(17) 
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Temperature 


FIG. 10: The l&2-cycle solutions of the Husimi lattice with J = — 1 and other parameters as 0. 


For a 2-cycle solution xi and X 2 we have 


Q+(x2,Xi) 

Q(X2,Xi) 

Q+(Xi,X2) 

Q(Xi,X2) 


(18) 


That is, the two solutions alternatively appear level by level. 

In this way, for a particular recursive relation, we can start from a set of initial seeds xi and X 2 to calculate the 
solutions on lower levels until we reach the recursively repeating solutions, which is called fix-point solution. Usually 
many initial seeds are tried to obtain all the possible solutions for a particular recursive relation. According to our 
experience, in the lattices discussed in this work, solutions with v > 3 were never obtained, while solutions with v = 1 
and V = 2 (1-cycle and 2-cycle solutions) are almost always available. Based on the property of x, which determines 
the probability that one site is occupied by the plus spin, a cycle solutions with w > 3 is hard to imagine. An example 
l&2-cycle solutions of the Husimi lattice, with J = — 1 and other parameters as 0 on a wide temperature range is 
shown in Fig. (TUI 

In Fig. [TUI at high temperature both 1-cycle and 2-cycle solutions are 0.5 and every site in the lattice has a 50% 
probability to be occupied by a -I- spin. This obviously represents a disordered state. At low temperature, the 1-cycle 
solution is still 0.5, while the 2-cycle solution occurs a transition and gives two branches, one of which goes to 1 and 
the other goes to 0 with temperature decrease. This indicates, for 2-cycle solution, that two neighbor sites will prefer 
to have different spin states. In the region close to zero temperature, if one site has 100% {x = 1) probability to be 
occupied by -I- spin, then its neighbor sites will have 0% {x = 0) probability to be occupied by -I- spin (i.e. 100% 
probability to be occupied by — spin). Therefore, at zero temperature we will have a plus and minus spins alternative 
arrangement on the lattice. Recalling the anti-ferromagnetic model we employed, this 2-cycle arrangement has the 
lowest energy (the most stable state) and corresponds to the ordered state (crystal). While the 1-cycle solution at 
the same temperature refers to the metastable state, that is, it is still a stable solution however with a higher energy. 
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The temperature where the 2-cycle solution appears is called the critical temperature, or the melting temperature 
Tni. At this temperature the amorphous state turns to be the ordered state (crystallization) or it can continue as 
a metastable state, the analog of three states here with liquid, crystal and supercooled liquid implies that it is the 
melting transition. The thermodynamic details will be discussed in section III. 


C. Recursive calculations on the surface lattice 

1. The calculation in the bulk 

For an infinite Husimi lattice, the calculation to get the fix-point solution is straightforward [25, 37], we can simply 
start with two artificial initial seeds and recursively calculate the solutions by Eq. [THImany times until we reach the 
fix-point solution. While for the SRL we developed in this paper, the bulk tree is finite and confined within surface 
trees, we need to carefully monitor the calculation on each specific site. 

Let us take a SRL with thickness of 2m + 3, we label the surface square as S, then the square next to S is labeled 
as the TO-th level, then origin square is labeled as 0. Firstly we start from two initial guesses xT X 2 on the surface 
sites as shown in Fig. [TTl The calculation through the surface square S provides the solution Xm+i on the top site of 
level m-th square. On the other two surface squares labeled as St and S” we use the same initial guess seeds again, 
however rotate their positions by putting aq on the top site and X 2 on two side sites, then the calculation based 
on S/ or 5"’ will give us the solution x/m+i and a:”m+i on the side sites of level m-th square. This is because for 
the anti-ferromagnetic Ising model and the properties of 2-cycle solution introduced above, we expect a 2-cycle style 
solution on our lattice (and the 1-cycle solution is just a special case when the two 2-cycle solutions are the same), 
thus the rotation avoids us to get the same results on neighbor sites, which makes the 2-cycle solution impossible. 

Consequently, with the solutions Xm +2 and Xm+i we obtain Xm on the top site of square level m — 1, and on the 
other branch the rotation of Xm +2 and Xm+i will give us Xm-i, then we can continue the calculation to the origin 
square 0 as shown in Fig. HHa). 

Once we reach the solution Xi and xti at the origin square, we have the states inside the bulk if the solution is the 
fix-point one. This part of calculation is called downward calculation, which is from the surface to the bulk origin. 
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Downward Calculation 
(from surface to bulk) 


Upward Calculation 
(back to surface) 


FIG. 12: The calculation direction though the bulk tree: (a) start from surface to the bulk origin, (b) start from the bulk origin 
and go back to surface. 


However, this set of bulk solutions is from the initial guesses on the surface; it may not be the solution of the real 
state. We still have to find the fix-point solution on the surface, and the bulk calculation from the surface fix-point 
solution is then the final result we are looking for. To find the fix-point on the surface, we start from the bulk origin 
and trace back to the surface (upward calculation), which provides a bulk contribution Tb to the recursive calculation 
on the surface. 

The upward calculation is shown in Fig. I12f bb At the origin square, from the contributions of three identical 
trees and the local weight of the origin square, we can calculate the solution on the lower site, which is labeled as 
= /(a;l,a;/l,a:”l,w;(r)). (Index U stands for ’upward’). Then from , xi2 and x"2 we can get x)^ and so on. The 
upward calculation will provide the solution x})^_^_l at the base site of the surface square unit. This solution will be 
used as x-q representing the bulk effect to surface in the surface calculation. 


2. The calculation along the surface 

The situation on surface is more complicated than in the bulk. From the previous introduction on downward 
and upward calculation, we can see a hint that depending on the direction of calculation, the solutions might be 
different on the same site, for example the Xm and x^ are different even they are on symmetric sites to the origin. 
This direction issue does not affect us to explore the solutions describing the bulk, however it is critical in surface 
calculation. Therefore, we firstly classify two directions on the surface with specific labeling, as shown in Fig. 1131 
In Fig. [T31 the labeling on the surface is in this way: the solution on the base site is xb, which is the from 

upward calculation; the solution on the sides of the square is label as aq, and on the top is x^] the prime marks the 
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FIG. 13: Two calculation directions on the surface of SRL. 



direction that goes out from the square, while the label without prime is the direction going into the square. The 
scheme A shows a “side-to-top” direction: the solution going out from the top site, T 2 ', is calculated from the solution 
going into the side site TT and the base site solution xb ■ Then from we go through a single bond to get the next 
xT. The scheme B shows a “side&top-to-side” direction: the solution going out from the side site, xl', is calculated 
from the solution going into the side site xi, the solution going into the top site X 2 , and the base site solution xb. 
Again from xi' we then go through a single bond to get the next X 2 - Only after we have done the calculations in 
both directions, we can get the surface solutions pair, the solution going into the top site X 2 , and the solution going 
into the side site xT, to do the subsequent bulk calculation. (Note that initially this pair is a guess we made to do the 
first iteration’s bulk calculation). The surface calculation scheme is shown in Fig. [TT] 

The upward bulk calculation provides us the bulk tree contribution and the solution xb on the base site of surface 
square. This solution can be taken as a constant in the surface recursive calculation. In Fig. 1141 we start from the 
Xb and the initial seed xT to do two calculations as mentioned. The recursive calculation in process A will eventually 
give us a fixed xT, while the process B will give us a fixed X 2 . This new fixed set of xT and X 2 is then the seeds we 
are going to use for the next iteration’s bulk calculation. Note one difference in process A and B is that in process 
A we only need xT for calculation, which is updated step by step, however in process B after each step we will only 
have a updated xT', while xT is a constant in calculation. Thus, in practice we always do the process A first to get a 
fixed xT, then use this xT as a constant in the calculation of process B. 

The calculation on the square unit, for example, xT^ = /(aq, X 2 , xb) in process B, is the same as the bulk calculation. 
The calculation through the single surface bond, for example, X 2 = /(aq') in process B, is also similar with the square. 
The difference is that the PPF only has two terms since there are 4 configurations of a single bond structure: 


Zmi + ) = ^^Zm+l{Sm+l)w{T), 
r=i 
2 

Zm{-) = y^^m+l(»S'm+l)w(r), 


r=3 


where the sum is over F 
polynomials: 


= 1 and 2 for Sm = +1 , and over F = 3 and 4 for Sm = • Then we have the 
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FIG. 14: Surface calculation scheme along the surface of SRL. 


2 

^ ^ (5*771+1 )'^(r), 

r^i 
2 

Qm—{ym-\-l) — ^ ^ ■^m+l(577i-^i)'Uj(r), 
r-3 

Qm(3^m+l) — Q?7i+(^m+l) “t“ Qm—(ym+l)- 
The ratio is the function of the ratio on higher level: 


Qm{Xm+l) 

The calculation scheme to reach the fix-point solutions both on the surface and in the bulk is shown in Fig. [THl We 
do an embedded recursive calculations to get the final solutions on the surface and in the bulk. The n iterations give 
us recursively updated solutions in the bulk until we find the fix-point solution, and a recursive calculation along the 
surface is embedded in one iteration to take the effect of updated bulk contributions each time. This process is to 
make sure that we counter the mutual effects from surface to bulk and from bulk to surface. Another necessity of this 
complicated process is that the bulk tree has a finite size thus the fix-point solution is not guaranteed to be reached 
during one downward calculation. For the structure with thickness > 19, we can always obtain the fix-point solution 
at bulk origin, although on the first several layers close to surface we can only have numerical calculations instead of 
exact solutions. (This ‘numerical depth’ depends on the energy parameters setting and it shows always being smaller 
than 19 according to our experience.) 
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fixed 
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nth iterations' bulk and surface 
fix-point solutions to be the 
same ("fix fix-point solution") 


FIG. 15: The calculation process scheme to reach the fix-point solutions in SRL. The vertical process is calculations through 
bulk tree and the horizontal process is along the surface. 


III. THERMODYNAMIC CALCULATION 
A. General free energy calculation method: the Gujrati trick 

Since our lattices are infinitely large, it makes no sense to calculate the partition function or free energy for the entire 
lattice. However an exact treatment called Gujrati trick has been well developed to deal with the thermodynamic 
calculation on recursive lattices in our group’s previous work [ 2 ^ [s^, IH, . By this technique we can approach the 
thermal properties in a local area (per site) by the partial partition functions Z{S) and solutions x we discussed in 
the previous section. 

The Helmholtz free energy is a function of the temperature and partition function Z\ 


F = -TlogZ 

Here we take a regular Husimi lattice as an example to introduce the thermal calculation derived by figuring out 
the partition functions. Recall Eq. 1101 since Zq counts the contribution of the whole system, it must be the sum of 
contributions of the six sub-trees on level 1 and 2, and the weight of the two local square units on the origin. If we cut 
off the two sub-branches on level 2 and joint them together, we will have an identical but smaller lattice. Similarly 
the partition function of this smaller lattice is: 


^2 = Zl{+) exp(/3i7) -b Zl{-) exp(-/3i7). 


We also cut off the 4 sub-branches on level 1 and form two smaller lattices and their partition function will be 


= Zi(-b)exp(/3H) -b Zl{-)e-xp{-j3H). 


As an extensive quantity the free energy of the whole system is the sum of the free energies of these three smaller 
lattices and the two local squares: 
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F{Zo) = 2F{Zi) + F{Z 2 ) + F(Ziocai) 


This yields 


T^local — ^^08(727 )■ 

Zi Z 2 

as the free energy of the two local sqaures, i.e. four sites (three paris of half-sites and the origin site). 
The free energy per site is: 




By substituting 




m 


and 


(1 ^m)i 


we have 


F = --riog(- 


+ {1 - XQ)'^e 


Bf\x\e^^ + (1 ~ + (1 ~ a:2)^e“^^] 


(19) 
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FIG. 17: The thermodynamic behavior of a Husimi square lattice with J = —1 and all other parameters to be 0. 


recall that for either 1-cycle or 2-cycle fix-point solutions we have Xo = X 2 , and 


Qo 


Bp 

BlB2' 


It follows 


F = --T\og{Q. 


^x\e^^ + (1 ~ xi)‘^e~^^^' 


With the calculation of Qp and fix-point solution xp and xi we discussed in the previous section, the free energy 
can be easily achieved. 

For a recursive structure with an origin point, we can always do this “cut and rearrange” trick to obtain the partition 
function and free energy per site around the origin point. 

The entropy is the first derivative of the free energy with respect to the temperature: 


S = -dF/dT. 

With the free energy and entropy we have the energy per site (energy density) as 


( 20 ) 


E = FFTS. (21) 

A typical thermodynamic behavior of a Husimi square lattice with J = — 1 and all other parameters to be 0 is 
shown in Fig. [T71 

In Fig. [13 the free energy of two solutions are the same at high temperature. At the melting temperature Tm 
the 1-cycle’s free energy differs from 2-cycle’s and is less stable. The continuous entropy at of 1-cycle implies the 
supercooled liquid state, while the entropy decrease of 2-cycle state implies the crystallization. Although here the 
entropy decrease is not a discontinuous jump as the typical behavior of first order melting transition, we may still 
treat it as the phases-coexisting point since the 2-cycle state is the most stable state here. 
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FIG. 18: Two example failure cuttings in SRL. 


The thermodynamics of 1 and 2-cycle solutions agree with our anti-ferromagnetic model. The anti-ferromagnetic 
Ising model prefer to anti-aligned neighbor spins, therefore the 2-cycle state should have the lowest energy, i.e. the 
crystal, while the 1-cycle solution represents the metastable state with higher energy. With the continuing decrease 
of temperature we can observe that at T = 1.1 the entropy of 1-cycle state (supercooled liquid) quickly decreases to 
zero and becomes negative. This is the Kauzmann’s paradox and the ideal glass transition. The negative entropy is 
unphysical thus the metastable state must undergo a transition to be the glass state at Tk, the ideal glass transition 
temperature. The free energy of the glass state below glass transition temperature is indicated by the last legend in 
Fig. [TT] Note that this branch is not provided by the calculation, instead it is a theoretical expectation. We can 
modify our calculation to achieve a stable solution of this glass state branch, but since this part is not our interest we 
have not done that calculation in this paper. 


B. The Gujrati trick applied on the SRL 

We are following the Gujrati trick to solve the free energy and consequent thermal functions of the SRL. Although 
the basic principle is the same, the asymmetrical structure of the surface lattice requires further complex tricks to do 
the calculation. If we take a random site on the surface as the origin, to do the “cut and re-joint” trick we need to 
select a local area around the origin, and find matching sub-trees contributing to the local area. Here the matchings 
are more specific than they are in a homogeneous bulk lattice. For example in the Husimi lattice discussed in previous 
section, all the sites are identical, thus if two sites have the same solution x on them, the sub-trees cut off from them 
can be joined together to make a new lattice. However, the sites on the surface of SRL have three different situations: 
on the top of a surface square where the sub-tree’s contribution going out, on the side of a surface square where the 
sub-tree’s contribution coming into, and on the bottom of a surface square where the bulk tree’s contribution coming 
into. The matchings must be done on sites with exactly the same situations to make an identical but small lattice, 
but the asymmetrical structure of SRL makes the cutting and matching impossible. Wherever we select the origin 
local area and cut the sub-trees, there will always be two sub-trees left and cannot be matched with each another. 

Fig. [TH] shows two failure cutting schemes, with “1 square and 2 single bonds” area and “2 squares and 4 single 
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FIG. 19: The modification on the surface of lattice and the cutting area around the origin. 


bonds” area respectively. In the cutting scheme ^1, by removing the local area surrounded by the discontinuous 
curves, we have four sub-tree at the point A, A’, B and B’. The two sub-trees on A and A’ can be linked together 
to form a smaller lattice with exactly the same structure. While the sub-trees on B and B’ do not match with each 
other, since they are a bulk sub-tree on B and a surface sub-tree contributing the top site of a surface square. 

In the cutting scheme #2, the local area of “2 squares and 4 single bonds” is bounded by the solid curves. The 
sub-trees on C - C’ and D-D’ can form smaller lattices. Although the new formation of D - D’ is not identical to the 
entire lattice, we can still count its partition function and handle the ratio calculation. However the same difficulty, 
as we encountered in scheme #1, presents on the sites E and E’. The E site is at the side of a surface square while 
the E’ site is on the top. 

Therefore, either we choose odd or even numbers of basic units as the local area around origin, the Gujrati trick 
cannot be done. We have to somehow modify the origin structure to avoid the problem of asymmetry. As shown in 
Fig. [T3l the single bond and square unit have a “head-to-side” connection, that is, following the calculation direction, 
the sub-tree coming from a single bond is always linked to the side site of next level’s square, or vice versa. This 
arrangement is designed to make the structure uniform on the surface, but it also causes the asymmetry. It is necessary 
to make one “head-to-head” connection as the origin of the surface, as shown in Fig. [THl 

In this way, we have two infinite and identical surface trees connected in a “head-to-head” style, this single bond 
is then the origin of our entire lattice, and except for this bond, all other surface single bonds are in “head-to-side” 
connection. Starting from the origin, the closest squares can be labeled as level 1; the second closest squares are on 
level 2 and so on. As shown in Fig. [T21 by the selection of four squares and six single bonds as the local area, we can 
form one identical lattice on level 1 at the sites A - A’, two identical lattices on level two at the pairs B - B’ and C - 
C’. The sub-trees cut at D, D’, E and E’ are identical bulk trees so we can pair them to form two bulk lattices. Now 
we can easily derive the free energy of the local area to be: 


Elocai = -7’log( 




( 22 ) 


Unlike the point origin in Husimi lattice, here the origin is a single bond unit. The total partition function thus 
has four terms due to the four possible states of the single bond: 
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FIG. 20: The redefined unit in SRL for the thermodynamic calculations around the origin. 


Zo = ZoU(+)zL(-)exp(-/3J) + ZoU(-)ZoL(+)exp(-/3J) + ZoU(+)Zo^(+)exp[/3(J + 27?)] + ZoU(-)zL(-)exp[/3(J-2H)] 

(23) 

To specifically track the spins in the origin area and rematch the sub-trees we define the two sub-trees meeting 
at the origin bond as “upper” and “lower” parts. In equation [521 the superscript U or L is to indicate the upper 
or lower part of the lattice. Similar integration can be employed to form the partition functions of smaller lattice 
we rearranged in Eal22l The ratio of partition functions can be represented by the solutions and polynomials which 
have been achieved in Section II. However in Section II, we go along the surface and calculate the solution on each 
site and the polynomial ratios step by step, these detailed solutions are unnecessary here and make it complicated to 
determine the Ea ]22l Hence we redefine the basic unit as one square and two single bonds linked on its side as shown 
in Fig. |20l 

By the selection of this rectangle basic unit, in the upper or lower branch we simply have two units contributing to 
the sides of next unit, and recursively the whole branch contributing to the origin bond, as indicated by the arrows in 
the lower part in Fig. [521 The contribution of bulk tree can be treated as a constant. The calculation based on this 
selection will provide only one solution as or x^ on the output site, regardless of it is 1-cycle and 2-cycle solutions 
on the surface. This single solution is sufficient for us to determine the Eg 1551 For instance, if we look at the smaller 
lattices rearranged on points A - A’ in Fig. [501 we will only need the solution x^ and x^ because they are on the 
output sites A and A’. 

Similarly, for either x^ or a:^, recall the Eg 1151 and Zm{+) = BmXm, Zm{—) = BmUm with Bm = [Zjn{+) -|-2’m(—)], 
then we have 




Bm+lBB 


2 

r {m+l} 


where the S'b is the index of the site linked to the bulk tree. 

With six spins the rectangle unit has 64 possible configurations thus we can introduce the polynomials: 
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Etn ('S'm+l )] (*S*B)^(r), 

F—1 {m+1} 
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Qm—(ym+l) — Etn 2m+l(5'm+l)]2B(5'B)w(r), 

F—33 {m+1} 


Qm(^m+l) — Qm+(^m+l) “t" Qm—iUm+l) — ^9 D 

-°m+l^B 

Then we can obtain the 1-cycle recursive relation: 


— 


Qm+ (^m+1) 
Qm{^m+1 ) 


And the fix-point solution is the one cycle solution or x^, Although the actual solution site by site on the 
surface could be either 1-cycle or 2-cycle. It is important to clarify that the solution x^ or x^ are exactly the same 
to the solutions we achieved on the surface in Section II. The reason we re-select the basic unit is to obtain a simple 
polynomial Q for the calculation in equation (3.13) and (3.14). 

Now we may rewrite Eql^^lwith H = 0 


^locai = -Tlog[ 




{x^y^ exp(—/3 J) -|- y^x^ exp(—/3 J) -I- x^x^ exp pj + y^y^ exp /3 J)^ 

1 


(24) 


(x^a;}^exp(/3iJ) -hi/^yj^exp(-/3i/))(a:}^a:]^ exp(/3i?) -H exp(-/3H))^ 
This local area contains 14 sites, thus the averaged free energy on each site is 


^ ^^Alocal- 


An reference free energy behaviors of 1-cycle and 2-cycle solution of SRL with thickness =19, neighbor interaction 
J = —1, surface neighbor interaction J = —1 and all other parameters to be zero is shown in Fig. 

An interesting phenomenon can be observed from Fig. 3.6. The free energy of 2-cycle solution (crystal state) is 
higher than the free energy of metastable state between T = 1.65 to 2.2. This is different from the bulk behavior in 
Fig. [51] The cross point of 1 and 2-cycle’s free energy at lower temperature can be determined to be the melting 
transition. Above the melting temperature, our results indicate that the anti-aligned spins arrangement is less stable 
than the 0.5 solution. This behavior has not been investigated, yet it is not clear whether it is simply unphysical or 
it implies a “super-heated crystal” state. 

Below the melting temperature the behavior is similar to the bulk system. The 2-cycle solution represents the 
crystal state with a lower free energy and continues to decrease to the free energy of ideal crystal. While the free 
energy of 1-cycle solution, the metastable state, will reach a minimum point than bind back, this unphysical behavior 
implies the ideal glass transition. 

With equations dOl and [21] we can easily calculate the energy density and entropy from the free energy. The Fig. 
1211 shows the entropy derived from the free energy in Fig. HU In Fig. HD the black arrow on the entropy show the 
melting transition at the cross point of free energies (T = 1.65), where the entropy of 1-cycle solution will step to 
2-cycle solution’s entropy as a first order transition. The ideal glass transition temperature Tk can be clearly observed 
by the negative entropy of 1-cycle solution. The detailed discussion on surface thermodynamics and the surface effect 
comparing with bulk system will be presented in the next section. 
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FIG. 21: The free energy behaviors of 1-cycle and 2-cycle solution of SRL with thickness = 19, J = —1, J = —1 and all other 
parameters to be zero 



FIG. 22: The thermodynamic behaviors of SRL with J = —1, J = —1 and all other parameters to be zero. 


IV. RESULTS AND DISCUSSION 
A. Discussion on the solutions 

As introduced in previous chapters, the thermodynamic calculations are mainly based on the solutions x, i.e. the 
ratio of partial partition functions on SRL. The 1 and 2-cycle solutions of a simple anti-magnetic field Husimi case 
is presented in Fig. [101 and we have discussed how the melting transition, crystal state and metastable state can be 
indicated from the solutions. The reference 2-cycle solutions of SRL with J and J = —1, other parameters as 0, and 
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FIG. 23: The reference 2-cycle solutions of SRL with J and J =-l, other parameters as 0, and thickness = 19. 


thickness = 19 are shown in Fig. [531 

The solutions of bulk Husimi lattice from Fig. [T0]are also included in Fig. j^for comparison. With the temperature 
increase, the solutions on surface converge to 0.5 at lower temperature than the Husimi bulk solutions. This indicates 
a lower melting transition temperature on the surface, which is easy to understand with smaller coordination number 
and less interactions on the surface, the spins on the surface are easier to be anti-aligned with less energy. Also, unlike 
the symmetrical solutions we usually have, the 2-cycle solution on the surface is asymmetric. This is because of the 
hybrid structure on the surface, which defines two different circumstances on the surface square unit and single bond 
unit respectively. Naturally the spins on the surface square unit has the solution closer to the bulk solution, while 
the spins on the single bond unit has a less stable (closer to 0.5) solution due to the lower dimension. Depending on 
the initial seeds adopted for surface calculation, we may also have the other set of surface solutions symmetric to the 
one in [231 

The 2-cycle bulk solutions and are the solutions at the bulk origin (Fig. [T3|) . Since the Husimi 

lattice plays the bulk portion in SRL, we can expect the bulk origin solution to be identical to the solutions of Husimi 
lattice if the thickness is large enough to ignore the effects of surface to bulk, and this expectation is confirmed in 
the region T < 2. However the 2-cycle bulk solutions differ from the Husimi lattice solutions, and converge to 0.5 
solution quickly above T = 2. Since the bulk solution comes to be almost steady with thickness larger than 19, this 
difference is not really caused by the surface effect. Our calculation requires a set of initial seeds for the recursive 
calculation as the procedure described in Fig. [T31 the bulk calculation takes the surface solution as its initial seed, in 
this way, once the surface solution reaches 0.5, the initial seeds of 0.5 will immediately affect the recursive calculation 
inside the bulk and converges the bulk origin solutions to be 0.5. The bulk origin solutions at converging point with 
different thickness are shown in Fig. 1241 We can see the convergence occurs with very slight differences no matter 
how large the thickness is, while theoretically we should have the bulk origin solutions to be identical with Husimi 
bulk solutions. Simple to say, because of the property of recursive calculation, the bulk origin solution will be lead 
away from the exact description by the effect of surface solution in the temperature region between surface melting 
and Husimi bulk melting transition temperatures, which is not really the “surface effect” in nature. 

This is also a reason we are not interested in calculating the thermodynamics of the whole SRL. Theoretically, the 
solution on each site of the SRL can be approached and therefore we can calculate the thermal properties. Nevertheless, 
there are three facts make this calculation unreliable. Firstly, as mentioned above, in the temperature region between 
surface melting and Husimi bulk melting transition, the bulk origin solution is affected by the seeds of 0.5 solutions 
on the surface. Secondly, the recursive calculation technique requires several steps to reach the fix-point solutions. 
This implies the calculations on the first few layers closing to the surface are numerical instead of exact calculation. 
Although for large thickness bulk this error can be neglected, the thermal behaviors associated with exact different 
layers is not useful. Thirdly, the recursive structure of SRL proportionally generates more surfaces with increasing 
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FIG. 24: The bulk origin solutions (upper branch) at converging point with different thickness. The curves of thickness = 19 
and 33 almost overlap and cannot be distinguished in the graph. 


thickness. Unlike the regular lattice, in which the contribution of surface will be neglected with a sufficient large 
bulk, the SRL will have half sites on the surface with infinitely large bulk trees. In this way, the thermodynamics on 
each site is just the averaged value of surface and bulk values. For a short summary, our approach of SRL is good to 
track the thermal behaviors on the surface with the account of bulk contributions, but not to discover the thermal 
behaviors change with various thicknesses. 


B. The transition temperatures reduction 

In section I we have reviewed the findings that the presence of a free surface dramatically decreases the transition 
temperatures of bulk system. By comparing the thermodynamics on the surface/thin film we achieved in section II 
and III and in the bulk system, our calculations also clearly indicate the reduction of both melting and ideal glass 
transitions temperature on the surface. Fig. [25] show the free energy comparison of Husimi bulk system and SRL: 

In SRL the melting and ideal glass transition temperature are dramatically decreased comparing to Husimi bulk. 
In Section I we introduced the empirical equation to describe the temperature reduction with the change of the 
thickness. Since the thickness dependence of transition temperature reduction is not available in our methods, and 
our calculations focus on the transitions right on the surface, we may compare our results to the ratio of glass transition 
temperatures of bulk and the thinnest free-standing film in others’ works, either experimental or simulation results. 

The Tk reduction ratio of SRL is 0.89/1.1 = 0.809. In Forrest and co-workers’ work, the Tg of the thinnest PS 
film they made is 300K and the Tg of bulk PS is 369K, the reduction ratio is 0.813 [10]. In Torres and co-workers’ 
MD simulation, this reduction ratio of a free standing film is 0.24/0.3 = 0.8 [12]. In de Pablo and co-workers’ MC 
simulation, this ratio is 0.85/1.08 = 0.79 [23]. 

Because our results are from the reference case, which is very simple without any particular artificial parameters 
setup, the fact that our ratios are close to others’ results may only verify the validity or practice of our method. To 
particularly describe a real system, the setup of energy parameters is the most critical issue. The effects of energy 
parameters in our model will be discussed in later section. The fact that similar reduction can be observed in our 
small moleculars model implies that the lower transition temperature on surface/thin film basically originates from 
the dimension reduction and less interaction constraints. The chain-structure of polymer system may not be the main 
reason for transition temperature reduction, although it does play important roles to affect the transition properties. 
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C. The effects of energy parameters 

In this section we are going to study the effects of different energy parameters in SRL by the variation of one 
parameter with other parameters fixed. The thickness is always set to be 19 for providing a sufficient tree size and 
relatively shorter calculation time consumption. 

In Eq[6]we specified the diagonal interaction energy parameter on the surface to be Jp, differed from the diagonal 
interaction Jp inside the bulk. And similarly we have J , (J/) and J”. This specification is to enable us setup different 
interaction circumstance on the surface and provide more versatility in our model. However the effects of either Jp 
or its counterpart on the surface Jp are the same. Thus in this section we are going to simplify the case and setup 
Jp and Jp , J”and J” to be the same. Nevertheless, a variation with fixed J will be discussed in details, because the 
nearest-neighbor interaction J and J have a much larger weight in the Hamiltonian and play more critical roles in 
the system. The difference setup of J and J provides the simulation of surface tension which is critical to the surface 
properties. 


1. The diagonal interaction Jp 

In the Hamiltonian Eql^l the nearest-neighbor interaction J is negative by the definition of anti-ferromagnetic Ising 
model, the system prefers an anti-aligned spins arrangement to obtain a lower energy with a negative first term, 
then we can also observe that for an anti-aligned spins arrangement, the diagonal spins pairs are in the same state. 
Therefore, unlike the nearest-neighbor interaction J, a positive diagonal interaction will make the second term to be 
negative, which encourages the alternating arrangement and increases the transition temperatures (i.e. the crystal is 
more stable to melt). On the other hand, a negative Jp will compete with J and trend to destroy the ordered state, 
which decreases the transition temperature. Because we have four nearest-neighbor and two diagonal interactions in a 
square unit, the nearest-neighbor interaction outweigh the diagonal interaction. Also, there is no diagonal interaction 
in the surface single bond unit, which also lowers the contribution of Jp. The variation of Jp will only moderately 
change the transition temperatures on the surface; the overall thermal behaviors are similar to the reference case. 
The thermal behaviors with four different Jp values ±0.2 and ±0.4 are calculated and the melting and ideal glass 
transition change with Jp variation is summarized in Eig. [55] and table 2. It is obvious that positive Jp will make the 
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FIG. 26: The transition temperature variations with different Jp. 


system more stable and increase the transition temperatures, or vice versa. We also found that Jp cannot be larger 
than ±0.5 otherwise stable solutions can not be reached. 

Table 2. The transition temperature variations with different Jp 



Tn, 

Tk 


0.4 

1.95 

1 

1.950 

0.2 

1.82 

0.95 

1.916 

0 

1.70 

0.89 

1.910 

-0.2 

1.55 

0.88 

1.761 

-0.4 

1.38 

0.80 

1.725 


2. The quadruplet interaction J” 

The quadruplet interaction J” is a complicated term. By simply analyzing the fourth term in Hamiltonian Eq. (2.1), 
it is difficult to give out a clear expectation on the effects of J” since both system preferred or defective structures 
can give out either positive or negative values of fours spins product, then consequently either positive or negative 
values of J” can against the crystallization. Transition temperatures of J” = ±0.2, ±0.4, ±0.6, ±0.8,and ±1.0 with 
other parameters fixed to be 0 are shown in Fig. [271 

From the graph above, we can see that the parameter J” has a similar effect to Jp: positive J” will make the 
system more stable and increase both transition temperatures, and vice versa. We also found that, unlike the limited 
value of Jp, J” can be assigned with large value such as ±1.0 without destroying the stable solutions. The transition 
temperatures variation with different J” are summarized in table 3. 
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FIG. 27: The transition temperature variations with different J”. 


Table 3. The transition temperature variations with different J” 

J” Tm Tk TJTk 
1.0 1.89 1.00 1.89 
0.8 1.88 0.98 1.92 
0.6 1.84 0.96 1.92 
0.4 1.78 0.94 1.89 
0.2 1.76 0.92 1.91 
0 1.70 0.89 1.91 

-0.2 1.64 0.86 1.91 
-0.4 1.56 0.80 1.95 
-0.6 1.46 0.76 1.92 
-0.8 1.36 0.69 1.97 
-1.0 1.26 0.62 2.03 


One important difference between the effects of Jp and J” can be observed in table 2 and 3: Although both 
parameters act similarly in changing the transition temperature, with the J” variation the ratio of Tin/Tk is relatively 
constant unless J” is assigned with extraordinary values like ±1.0, while with the Jp variation the ratio of T„i/Tk 
also changes dramatically. This difference could be useful in modifying the parameters setup to describe the real 
systems or experimental results. For a real system, factors affecting the thermodynamics can act in different ways, 
some of them may lift/reduce both transition temperatures but keep the supercooled liquid region constant, or some 
may change the window size of supercooled liquid state. 


3. The surface nearest-neighbor interaction J 

As mentioned above, for convenience we setup the energy parameters in the bulk and on the surface, such as Jp 
and Jp, to be the same. Nevertheless, our lattice has a fixed bond length (this length is unit 1 in Ising model) 
between neighbor sites on the surface and in the bulk, and we also know that the asymmetric interactions and the 
surface tension are the main reason of unique properties of the surface. Therefore with a fixed bong length, we have 
to somehow modify the surface interaction circumstance to mimic the surface tension or asymmetric interactions. 
This can only be made possible with different J and J setup. Thus we fixed the bulk nearest-neighbor interaction 
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FIG. 28: The transition temperature variations with different J 


parameter J to be —1, and observe how the variation of surface bond interaction J affect the thermodynamics on the 
surface. The transition temperatures change with variation of J = —0.5, —0.7, —0.9, —1.1 — 1.3 is summarized in Fig. 
[28] and table 4. 

From Fig. |5S|we can conclude that larger absolute value of negative J makes the system more stable and increases 
both melting and ideal glass transition temperatures. The variation of J revealed that the surface is more stable 
with large surface tension (—J > — J), or easier to undergo transitions than in the bulk with small surface tension 
(—J < — J). Positive J was also been tested however no stable solution can be reached. This is because the positive 
J will prefer ferromagnetic-aligned spins arrangement which only has 1-cycle solution available. On the other hand 
it is also hard to imagine a homogeneous system with particles attractive to each other in the bulk but repulsive on 
the surface. _ 

Table 4. The transition temperature variations with different J 


J 

Tm 

Tk 

T 

m/Tk 

-0.5 

1.12 

0.55 

2 

.04 

-0.7 

1.38 

0.70 

1 

.97 

-0.9 

1.60 

0.82 

1 

.95 

-1 

1.70 

0.89 

1 

.91 

-1.1 

1.80 

0.94 

1 

.91 

-1.3 

2.00 

1.05 

1 

.90 

-1.5 

2.18 

1.14 

1 

.91 


4- The triplet interaction J’ and magnetic field H 

It is easy to understand that we have 0.5 solutions at high temperature, and 2-cycle solutions symmetric to 0.5 
at low temperature, is because that the magnetic field H is 0, the system has no bias on either -|- or — spins on a 
particular site. Thus at high temperature one site has 50% probability to be occupied by either -|- or — spin. If a 
non-zero value of magnetic field H exists, we can expect the 1-cycle solution off the central 0.5 line or asymmetric 
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2-cycle solutions. A positive H will prefer more -I- spins and moves the 1-cycle solution higher than 0.5, or vice versa. 
In our previous research on bulk Husimi lattice, the three spins interaction (triplet) Jf acts similarly to magnetic 
field H in changing the bias of spins. However, in SRL the J/ and magnetic field H are not useful, although for a 
general presentation we still have them in our lattice model and energy equations. The reason is that, the uneven 
property of the surface on SRL determines the 1-cycle solution can only be 0.5. If a non-zero magnetic field H or 
triplet interaction J/ presents, the 1-cycle solution on the surface cannot be achieved. The single bond unit does not 
have triple interaction J/ and the effect of H is much smaller than it is in the square unit, therefore the anti-aligning 
property of the single will always prefer to have different spins on it. In another word, if a 1-cycle solution can be 
achieved on the single surface bond, it must be the 0.5 solution. Even we can still calculate the thermodynamics of 
2-cycle solution on the surface, we must have the corresponding 1-cycle solution to determine the melting and ideal 
glass transitions. Therefore, we have to abandon the variation of H and J/, which is a disadvantage of the SRL. 


V. CONCLUSION 

We applied recursive lattice technique to investigate the thermodynamics and ideal glass transition on the sur¬ 
face/thin film of Ising spin system with a theoretical base. A recursive lattice (SRL) were constructed to describe a 
regular finite size square lattice, with a 2D bulk surrounded by ID surface. The SRL utilizes the finite-sized Husimi 
lattice, which is integrated by 2D square units, to be the bulk part, and adopts single bond to connect the finite-sized 
bulk parts. These single bonds and the outsider bonds of bulk trees assemble the surfaces surrounding the bulk 
part and slipping through independent bulk parts. The lattice is constructed with particular approximations and 
compromises and is believed to be reliable approximation to the regular square lattice in the aspect of coordination 
numbers, i.e. the number of neighbor sites is 4 inside the bulk and 3 on the surface. 

The Ising spins of two possible states -|- or — are put on the sites of SRL. We assigned anti-ferromagnetic interaction 
in the model and different spin states in neighbor pairs are preferred to be the stable state (crystal). The partition 
function of a sub-tree with its base spin state fixed is defined as partial partition function (PPF). Based on the 
recursive properties, the PPF of one sub-tree can be expressed as a function of the PPFs of its sub-sub-trees and a 
local weight. This recursive relation enables us to derive the ratio of PFFs on one site. This ratio is called solution 
from which the thermodynamics of system can be achieved by Gujrati trick. Two kinds of stable solutions usually can 
be achieved by recursive calculation. One is in 2-cycle form and represents the ordered state, the other is in 1-cycle 
form representing the amorphous/metastable state. 

The thermal behaviors of 2-cycle and 1-cycle solutions indicate the melting transition by the cross point of free 
energies of two solutions, and the ideal glass transition by the negative entropy of 1-cycle solution. Our results agree 
with others’ work that the transition temperatures are dramatically reduced on the surface/thin film comparing which 
in the bulk. Nevertheless our work shows that, regardless of specific properties of particular materials, either polymer 
or small molecules, this transition temperature reduction can be simply due to the dimension downgrade and lower 
interaction energy on the surface. 

The effects of different interaction energy parameters other than nearest-neighbor interaction such as diagonal 
interaction and four spins (quadruplet) interaction were investigated. The variation of energy parameter could either 
increase or decrease the stability of system and change the transition temperatures according to the Hamiltonian. 
The behaviors of energy parameters can be used to setup a particular model to describe the real system. 

While our finding agrees well with others’ experimental and simulation results, comparing to others’ work, there are 
two advantages in our approach: 1) in this work our model focuses on the small molecules system, it reveals the basic 
dimension origin of transition temperatures reduction without involving the long chain properties of polymer system; 
2) the thermodynamics of systems are derived by exact calculation method, the computation time is much shorter 
than typical simulation methods, usually the calculation of one set of parameters in the interesting temperature region 
can be done in less than 100 seconds. 
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